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We numerically study the ground state of an S=l/2 antiferromagnetic Heisenberg model on a 
spatially anisotropic triangular lattice, which is an effective model of Mott insulators on a triangular 
layer of organic charge transfer salts or CS2CUCI4. We apply a numerical variational method by using 
a tensor network with entanglement renormalization, which improves the capability of describing a 
quantum state. We identify magnetic ground states for 0.7 < J2/J1 < 1 in the thermodynamic limit, 
where Ji and J2 denote the inner-chain and inter-chain coupling constants, respectively. Except for 
the isotropic case (Ji = J2), the magnetic structure is spiral with an incommensurate wave vector 
that is different from the classical one. The quantum fluctuation weakens the effective coupling 
between chains, but the magnetic order remains in the thermodynamic limit. In addition, the 
incommensurate wave number is in good agreement with that of the series expansion method. 



I. INTRODUCTION 

The physics of Mott insulators has been attracting at- 
tention since the discovery of high-temperature supercon- 
ductors. In the past decade, a number of new Mott insu- 
lator materials on a triangular layer have been found. For 
example, Cs2CuCl4^ and organic charge transfer salt^, 
such as /^-(BEDT-TTF)2 X and /3'-Z[Pd(dmit)2]2, have 
been extensively studied using experimental and theoret- 
ical approaches. At low temperatures in the Mott insula- 
tor phase, these materials show various quantum states: 
an antiferromagnetic long-range ordered state, a valence 
bond crystal state, and a disordered state. In partic- 
ular, the disordered behaviors in /^-(BEDT-TTF)2 Cu2 
(CNja'^, EtMe3Sb[Pd(dmit)2]J^, and Cs2CuCl4^ are 
of great interest. 

The simplest effective model of spin degrees in Mott 
insulators is the 5 = 1/2 antiferromagnetic Heisenberg 
model on a triangular lattice. Since the triangular layer 
in a real material is distorted, we have two groups of 
Heisenberg interactions^! In Fig. [ija), we show these 
two kinds of interactions as solid and dotted links. The 
Hamiltonian is written as 
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where (ij) and {{ij)) denote pairs of sites on solid and 
dotted links in Fig. [TJa), respectively. The coupling co- 
efficients Ji and J2 are positive. The ratio J2/J1 in 
real materials varies widely, from | to 1. For example, 
J2/J1 for /^-(BEDT-TTF)2 Cu2 (CN)3 was estimated to 
be close to 1, and that for CS2CUCI4 was estimated to be 
about 

The model of Eq. ([T]) interpolates among indepen- 
dent chains (J2 = 0), the fully frustrated triangular 
lattice (Ji = J2), and the unfrustrated square lattice 
(Ji = 0). Geometrical frustrations are present through- 
out the model, except at two special points (Ji = and 
J2 = 0). In the classical case, the ground state can be 
solved exactly: There are two long-range order phases 
at zero temperature [see Fig. [ijb)]: a Neel state on a 
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FIG. 1. (Color online) (a) Spatially anisotropic triangular 
lattice with two groups of interactions, denoted by dotted 
and solid links. Ji and J2 denote the coupling coefficients on 
each link, (b) The magnetic moment twist angle along the 
J2 axis of the classical model (solid line) and the proposed 
phase diagrams of the quantum model (horizontal strips) on 
the spatially anisotropic triangular lattice. Here, SP and SL 
denote the spiral and spin- liquid phase, respectively. 



square lattice for J2/J1 > 2 and a spiral state with a 
smoothly changing wave number for J2/J1 < 2. How- 
ever, the phase diagram in the quantum model cannot 
be solved analytically. Thus, we have used various nu- 
merical or approximate methods. First, the ground state 
at the isotropic point (J2/J1 = 1) has been studied. For 
example, both exact diagonalizatioiP and density matrix 
renormalization group (DMRG) calculations^ reveal a 
120° magnetic ordered ground state, whose wave number 
is equal to that in the classical model. Although quantum 
fluctuations reduce the magnetic moment magnitude, the 
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FIG. 2 
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(Color online) (a) Matrix product state for six sites: 
Tsi^.^-,SQ- (b) Tensor network with entanglement renormaliza- 
tion on a triangular lattice. Solid (red) lines represent tensor 
contractions for two connected tensor indices. 



120° magnetic ordered state was confirmed. However, the 
stability of the magnetic ordered state in the anisotropic 
region is controversial. Yunoki and Sorrella^^ reported a 
disordered state for J2/J1 ^ 0.79 by using a variational 
Monte Carlo (VMC) technique. In addition, Heidarian et 
al.^^ reported the disappearance of magnetic long-range 
order at J2/J1 ^ 0-85 by using another VMC technique. 
Weng et al)^ also reported a similar disordered state for 
J2/J1 < 0.78 by using the DMRG method. Reuther and 
Thomald^ reported a disordered state with collinear an- 
tiferromagnetic stripe fluctuations for J2/J1 < 0.7 ~ 0.9 
by using the pseudofermion functional renormalization 
group method. Thus, the disordered behavior in real ma- 
terials may be captured by these states. However, some 
reports are contradictory, as shown in Fig. [ijb) . Zheng 
et al. proposed a spiral phase for < J2/J1 < l-H by 
using a series expansion methocf^. Weichselbaum and 
Whit^IS also reported a long-rang magnetic correlation 
with an incommensurate wave vector in the whole region 
of < J2/J1 < 1 by using the DMRG method with dif- 
ferent boundary conditions. The renormalization group 
analyse d^^'^^ suggest a direct transition from spiral to 
collinear antiferromagnetic order at J2/J1 ^ 0.3. There- 
fore, the stability of the spiral state in the quantum case 
is crucial for understanding the physical behavior of real 
materials. 

We used a new numerical approach to calculate the 
ground state. Usually, quantum Monte Carlo (QMC) 
methods are powerful tools for two-dimensional quantum 



models, because they are unbiased. However, the weight 
of QMC samples can be negative in frustrated quantum 
magnets, leading to a cancellation in sign, and the ac- 
curacy of simulations fatally decreases. Exact diagonal- 
ization can only be applied to small systems. Thus, we 
have chosen a variational method. In particular, the key 
point of our calculation is the trial wave function. It is 
based on a tensor network with entanglement renormal- 
ization (ER)^^. A tensor network is a theoretical tool 
in the field of quantum information to describe a quan- 
tum state. By modifying the network structure, we can 
freely design the structure of entanglements that mean 
quantum correlations in a quantum state. In general, 
the entanglement entropy of a subsystem is proportional 
to the area of the boundar};^^. A tensor network with 
ER also obeys the area law of entanglement entropy. 
Though it only has a bias owing to the particular net- 
work structure used in the calculation, systematic error 
can be controlled, in principle, by increasing the dimen- 
sions of tensor indices. Therefore, the tensor network 
method is regarded as one of the most promising tech- 
niques for treating numerically unsolved problems such 
as the present one. Unfortunately, successful applica- 
tions to quantum frustrated magnets in two dimensions 
are very fev^^. In what follows, we demonstrate the 
usefulness of ER by applying it to the model of Eq. ([T]) 
to clarify the nature of its ground state. 

By using the ER tensor networks shown in Fig. |2|b) , 
we confirmed the spiral state with incommensurate wave 
numbers for 0.7 < J2/J1 < ^ that overlaps with those of 
the disordered (spin liquid) phase reported in previous 
workJ™^ [see Fig. [l](b)]. In our results, quantum fluc- 
tuations weaken the effective coupling between chains, 
but the long-range magnetic order remains in the ther- 
modynamic limit. In addition, the incommensurate wave 
number in our results is in good agreement with that ob- 
tained by the series expansion method^^. Since these 
two approaches are different, our results provide strong 
evidence for the stable spiral phase. 

The paper is organized as follows: In Sec. [ll| we will 
briefly introduce a tensor network with ER designed for 
triangular lattice models. In Sec. |III[ we will report 
numerical calculations of the S = 1/2 antiferromagnetic 
Heisenberg model on a spatially anisotropic triangular 
lattice. In Sec. 



IV, we will summarize our results. 



II. TENSOR NETWORK WITH ER ON A 
TRIANGULAR LATTICE 

A. Tensor network 

Formally, the probability amplitudes of a wave func- 
tion can be regarded by a rank- TV tensor T as 



where N denotes the number 



of sites. However, we cannot treat a large- TV-site sys- 
tem by using a tensor, because the number of elements 
in a rank-A^ tensor exponentially increases. To avoid the 
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FIG. 3. (Color online) Coarse-graining transformation for a 
triangular lattice model. 



exponential increase, we replace the original large-rank 
tensor by using a set of tensor contractions of small-rank 
tensors. 

The tensor contractions can be drawn as a network. 
Thus, this type of wave function is called a tensor net- 
work wave function or, simply, a tensor network. The 
node of the network denotes a tensor, and the leg of the 
node denotes a tensor index. An edge connecting two 
legs represents a tensor contraction for the two corre- 
sponding indices. For example. Fig. [2|a) shows the ten- 
sor network for the probability amplitude as T^^^^^^ = 
Y.ti,-,te ^tQ.s^M ' ' ' ^t5,s6,t6- The one-dimensional ten- 
sor network is called a matrix product state (MPS). It 
is the general form of the wave function used in DMRG 
calculations. 



B. ER on triangular lattices 

Various types of tensor networks have been proposed 
for many-body quantum systems. The structure of the 
network affects entanglements in a tensor network state. 
For example, the MPS breaks the area law of entangle- 
ment entropy in more than two dimensions. Thus, in 
principle, it is not suitable for capturing the quantum 
state in two-dimensional quantum systems. To construct 
a tensor network for a triangular lattice model, we use a 
coarse-graining transformation removing short-range en- 
tanglements between coarse-grained regions. This is the 
ER method proposed by Vidal^^. In particular, since the 
tensor network with ER obeys the area law of entangle- 
ment entropy, it can describe a quantum state with large 
entanglements in principle. 

No systematic studies have been conducted on the op- 
timal network structure. The empirical rule is that it 
should decrease the entanglement between coarse-grained 
regions as much as possible and, at the same time, keep 
the computational cost of tensor contractions reasonable. 
We found a suitable ER tensor network for triangular lat- 



tice models, as shown in Fig. |2jb). It transforms a trian- 
gular lattice [Fig.jsja)] to a coarse-grained one [Fig.jsj^f)]. 
After the transformation, the number of sites decreases 
by a factor of 19. The coarse-grained unit cell is the filled 
gray hexagon in Fig. [3|b). As shown in Fig. ^h), the 
network consists of three sublayers. Each sublayer is oc- 
cupied by a single type of tensor: ix(red), 'u(green), and 
iL'(blue) from the bottom sublayer to the top sublayer, 
respectively, as shown in Fig. [2|b) and Fig. [s] Tensors 
u and V are called disentanglers^ because their purpose 
is to decrease short-range entanglements between coarse- 
grained regions. The tensors have upper and lower legs. 
An upper leg in one sublayer is connected to the lower 
leg of a tensor in the higher sublayer. Disentangler u has 
three lower legs and three upper ones. Disentangler v has 
six lower legs and two upper ones. Tensor w transforms 
seven sites to one site; this process is called isometry. In 
principle, the ER can be applied iteratively, and we usu- 
ally finish the ERs corresponding to the top tensor on 
the last coarse-grained lattice, which is a simple isome- 
try. In particular, the tensor network with multiple-level 
ERs is called a multiscale entanglement renormalization 
ansatz (MERA). 



C. Computational costs of MERA 

Ah tensors in MERA are isometric: {TIYtI = 5ij, 
where Tj denotes the tensor's element with index i (j) 
of upper (lower) legs. Because of the isometric property, 
an expectation value of a local operator can be evaluated 
on a subnetwork that is finite and much smaller than the 
whole network, in most a pplications. This subnetwork is 
called a causal con^^^^. Figure [Z] shows a causal cone 
for the expectation value of an operator on a triangle pla- 
quette of nearest neighbor sites. The number of tensors 
in a causal cone is proportional only to the logarithm of 
system size. Thus, the computational cost depends on 
the system size only weakly, compared to the exponen- 
tial growth that is naturally expected. It only increases 
by a polynomial of dimensions of tensor indices. 

We assume that tensor legs at the same "height" have 
the same dimensions. Then the size of tensors in MERA 
can be specified by only an integer set as (xi7X2,X3) 
Fig. [2|b), where Xi is the dimension of the upper index 
of tensors in the ith sublayer. The computational cost 
of the expectation value of operators on a triangle pla- 
quette becomes a polynomial of these integers. As shown 
in Fig. [4j all tensors except the local operator drawn as 
(yellow) tensor A are paired in a causal cone. We first cal- 
culate tensor contractions between paired tensors. Then, 
the causal cone is transformed to a tensor network that 
has a half height. Each tensor in the new tensor network 
is defined by paired tensors in the original network. If 
and only if an edge connects unpaired tensors, the edge 
remains in the new tensor network, as shown by the solid 
(red) lines in Fig. [4j Thus, the shape is similar to the up- 
per part of the original one. The number of remaining 
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FIG. 4. (Color online) Causal cone for an expectation value of 
a local operator A on a triangle plaquette of nearest neighbor 
sites i, j, and k: {ip\A{i, j, k)\'ip) . Here \ip) is represented by 
the tensor network with an ER level defined in Fig.|2jb) and a 
top tensor that covers six sites as in Fig. [5] The local operator 
A is drawn as a (yellow) tensor that has three upper and 
three lower legs. The upper and lower parts of the (yellow) 
tensor are subnetworks from \ip) and {ip\, respectively. We 
only draw tensor contractions between unpaired tensors as 
solid (red) lines. A tensor's index without a solid (red) line 
in the upper part is always connected to that at the same 
position in the lower part. The width of a causal cone is 
defined by the number of solid (red) lines at the same height. 
The maximum width of this causal cone is six. 



lines at the same height in the new tensor network is 
called the width of the causal cone. Usually, we calcu- 
late this tensor network from the local operator drawn 
as (yellow) tensor A in Fig. |4j The maximum number 
of indices of intermediate tensors is roughly double the 
maximum width of the causal cone. Thus, the mem- 
ory size needed for calculating a causal cone rapidly in- 
creases with increasing dimensions of tensor indices. In 
fact, the maximum width of the causal cone in Fig. |4] 
is six. Since we use multithreaded subroutines for ten- 



sor contractions, the total memory size needed to cal- 
culate a causal cone is limited by the memory size of a 
computational node, which strictly limits the maximum 
dimensions of indices. In addition, the maximum polyno- 
mial degree for the computation of the causal cone is also 
larger than the double maximum width of causal cone. 
For example, the maximum polynomial degree in Fig. [4] 
is 14. Since the Hamiltonian of Eq. ([T]) is written as the 
summation of local Hamiltonians on triangle plaquettes 
of nearest neighbor sites, the main part of the variational 
method can be decomposed into calculations of indepen- 
dent causal cones corresponding to local Hamiltonians. 
Thus, this part can be perfectly parallelized. Our main 
calculation has been done using the facilities of the Su- 
percomputer Center, Institute for Solid State Physics, 
University of Tokyo. We used 256 nodes in the largest 
case in which the tensor network with two ER levels was 
used. 



III. NUMERICAL RESULTS OBTAINED USING 
A TENSOR NETWORK WITH ER 

A. Isotropic triangular lattice 

First, we calculate ground states of finite and infinite 
systems for Ji = J2. The 120° magnetic ordered state at 
the isotropic point has been confirmed by previous works 
(see Table III in Ref. |23]). The purpose of this calcula- 
tion is to test the variational wave function defined in the 
previous section and to see the behavior for an S=l/2 an- 
tiferromagnetic Heisenberg model on a triangular lattice. 



1. Tensor network 

The wave function consists of the single ER level in 
Fig. |3] with a top tensor. The top tensor covers six coarse- 
grained sites after the ER. Thus, this tensor network 
structure is applied to TV = 6 x 19 = 114 sites. Fig- 
ure [5] shows the tensor network structure. Large solid 
circles denote the positions of coarse-grained sites. We 
put the top tensor on the parallelogram frame in Fig. [5j 
We apply this tensor network structure to both finite and 
infinite lattices. In this paper, we call the former the pe- 
riodic boundary condition (PBC) scheme and the latter 
the infinite-size scheme, respectively. 

In the PBC scheme, the total number of sites is just 
6 X 19 = 114. We set a skew PBC so that ah parallelo- 
gram frames in Fig. [5] are the same. We notice that this 
PBC is consistent with a three-sublattice structure of a 
triangular lattice. In contrast, in the infinite-size scheme, 
we arrange the same tensor network structure, with 114 
sites per unit cell on an infinite lattice. Then the top 
layer is defined as the product state by top tensors. In 
addition, we assume that the tensors at the correspond- 
ing positions in all repeated units are the same. Thus, 
we can define a wave function by the finite set of tensors 




FIG. 5. (Color online) Tensor network structure with a single 
ER level and a top tensor for six coarse-grained sites. Big 
solid circles denote positions of sites on the coarse-grained 
lattice. The top tensor is put on the parallelogram frame. 
All parallelogram frames are the same by a skew periodic 
arrangement. 

for the infinite lattice. This type of MERA is called a 
finite-correlation MERA^^, because reduced correlations 
become exactly zero for large distances. The distance 
limit for finite reduced correlations is roughly the size 
of the unit cell. The main difference between the two 
schemes lies in the causal cones. In the PBC scheme, 
all the causal cones are limited to one unit cell, because 
of the PBC. However, in the infinite-size scheme, some 
causal cones extend into multiple unit cells. Thus, the 
computational time for the infinite-size scheme may be 
longer than that for the PBC scheme. 

We assume that all tensors are independent within the 
unit cell to have more variational freedom with less bias. 
We optimized them to minimize the total energy of the 
tensor network states. The tensors are iteratively up- 
dated by the singular value decomposition method^^. Al- 
though this minimization problem may have some local 
minimum states, we obtained stable results starting from 
random initial tensors. 



2. Energy 

Our ER tensor network has the spatial structure shown 
m Fig. |3l As we mentioned above, if we use the fi- 
nite dimension of the tensor index, the network-structure 
bias may cause a significant systematic error. To check 
whether or not this is the case, we try several sets of 
dimensions of tensor indices. Figure [6] shows local en- 
ergies on finite and infinite lattices for various tensor 
sizes. The value of the local energy on a triangle pla- 
quette defined by nearest neighbor sites i, j, and /c, 

• Sj + Sj • S/e + S/e • S^, is shown on a color scale. 
In the 120° state, which we believe is the ground state 
in the present case, the local energy is homogeneous. 
Therefore, we expect that the whole system should be 
uniformly colored if the error of the calculation is suffi- 
ciently small. Figures [6ja) and ^c) show the results of 
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FIG. 6. (color online) Local energies on the triangular plaque- 
ttes. The value of the local energy on a triangle plaquette de- 
fined by nearest neighbor sites i, j and k, S^-Sj+Sj -Sfc+Sfc-Si, 
is shown by the color scale. The bottom bar shows the cor- 
respondence between the value of local energy and color. Re- 
sults of the PBC scheme are shown in (a) and (b). Those of 
the infinite-size scheme are shown in (c) and (d). The size of 
tensors in (a) and (c) is (xi5X2,X3) — (2,2,2). In (b) and 
(d), (xi,X2,X3) = (2,8,8). 



small tensor size, (X17X27X3) = (2,2,2), for PBC and 
infinite-size schemes of TV = 114, respectively. There are 
clear spatial patterns resulting from the structure of ER 
in both cases. Increasing the tensor size should improve 
the quality of the tensor network states. In fact, as shown 
in Figs. |6|b) and[6]^d), the patterns are clearly more ho- 
mogeneous than for small tensors. We notice that the 
patterns in the infinite-size scheme are quite similar to 
those of the PBC scheme. Thus the assumption of direct 
product states for infinite systems does not affect the spa- 
tial pattern of local energies in the tensor network states. 

Figure [7| shows the energy per site, for the tensor 
sets (xi,X2,X3) = (2,2,2), (2,4,4), (2, 8, 4), and (2, 8, 8) 
for both PBC and infinite-size schemes. When the ten- 
sor size increases, the energy of tensor networks is in- 
deed improved. The lowest energy per site at (2,8,8) 
is EpBc{N = 114) = -0.54181 and Eir,f{N = 114) = 
—0.54086. These values compare well with results ob- 
tained from other methods (see Table III in Ref. |23|) . 
In particular, the result of a Greenes function quantum 
Monte Carlo (GFQMC) calculation with stochastic re- 
configuration (SR)24 is E{N = 144) = -0.5472(2). 
Direct comparison may be difficult, because our skew 
boundary is not equal to the periodic one used by the 
authors of Ref. 24 and our lattice {N = 114) is smaller 
than their lattice {N = 144). However, our result is 
close to their result (:1% lower than ours). Our result 
for the infinite-size scheme also agrees well with previ- 
ous estimates for the thermodynamic limit. In particu- 
lar, it compares favorably to that of a series expansiorP^, 
E = -0.5502(4), and to that from a GFQMC with SR 
calculation^^, —0.5458(1), in the thermodynamic limit. 
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FIG. 7. (Color online) Energy per site, and the average 
of entanglement entropies on triangle plaquettes, (Sijk), for 
PBC and infinite-size schemes (N — 114). Circles (red) and 
crosses (blue) denote E for PBC and infinite-size schemes, 
respectively. Triangles (green) and squares (yellow) denote 
{Sijk) for PBC and infinite-size schemes, respectively. The 
error bar of E shows the deviation of local energies on trian- 
gular plaquettes. 
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FIG. 9. (color online) Averages on-site magnetizations on 
triangular plaquettes. The bottom bar shows the correspon- 
dence between the value of average on-site magnetizations and 
color. The schemes and tensor sizes from (a) to (d) are equal 
to those from (a) to (d) in Fig. [g] respectively. 



the tensor size increases, the average value of entangle- 
ment entropies also increases, as shown in Fig. j?! How- 
ever, as shown in Figs.jsjb) andjsjd), the spatial patterns 
of entanglement entropies are different. In the infinite- 
size scheme, the boundary of unit cell has weak entangle- 
ment entropies. The reason for this is that the top layer 
is a direct product state. Although the assumption of a 
direct product state does not affect the spatial pattern 
of local energies, that of entanglement entropy is more 
sensitive. Thus, the entanglement entropy may be useful 
for checking wave function quality in other cases. 
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FIG. 8. (color online) Entanglement entropy of triangular 
plaquettes defined by nearest neighbor sites i, j, and k in ER 
tensor network states. The bottom bar shows the correspon- 
dence between the value of entanglement entropy Sijk and 
color. The schemes and tensor sizes from (a) to (d) are equal 
to those from (a) to (d) in Fig. p] respectively. 



3. Entanglement entropy 

Entanglement entropy measures the quantum corre- 
lation between a considered region and another region. 
Figure |8] shows the entanglement entropy of a triangular 
plaquette on a color scale. It is defined as 

Sijk = -T^^iPijk^^Pijk], (2) 

where pijk is the reduced density matrix of sites i, j, and 
on a triangle plaquette. There is clear spatial inho- 
mogeneity owing to the network- structure bias, as is also 
seen in the local energy. The entanglement entropies on 
the boundary of coarse-grained regions are lower. When 



The magnetization on a site i is defined as 

Mi = V(S,) ■ (Si), (3) 

where (•) denotes the expectation value of the operator 
by a variational wave function. Figure [9] shows the av- 
erage on-site magnetizations on a triangular plaquette. 
They depend on the tensor size, and they show spatial 
patterns that depend on the structure of the ER tensor 
network. In addition, as seen in the case of entanglement 
entropy, the direct-product nature on the top layer also 
gives rise to spatial inhomogeneity. To make the extrap- 
olation possible, we plot estimates of the magnetization 
as a function of the corresponding estimates of energy. 
Figure [lO] shows the average on-site magnetizations ob- 
tained by using the infinite-size scheme (A^ = 114), M, 
for various tensor sizes. The vertical and horizontal axes 
denote M and Einf{N = 114), respectively. The depen- 
dence of M on the tensor size is high. While the mag- 
netization for the largest tensor size is 0.322(2) [see the 
left-most (red) circle], we cannot simply extrapolate it. 
The solid (red) curve is fitted to points of M. If we 
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All values in Fig. 10 are 120.0(4). Therefore, the ground 
state of an S=l/2 antiferromagnetic Heisenberg model on 
an isotropic triangular lattice is a magnetic ordered state 
with 120° structure. This result is consistent with results 
from many previous works (see Table III in Ref. | 



FIG. 10. (Color online) Average on-site magnetizations and 
average angles between magnetic moments of nearest neigh- 
bor sites. These results are obtained by the infinite-size 
scheme (N = 114). The horizontal axis denotes energy per 
site. The left and right vertical axes denote the average on- 
site magnetizations on all sites and at the centers of ER, 
M and Mc, and the average angle between magnetic mo- 
ments of nearest neighbor sites, 0, respectively. Circles (red), 
squares (green), and triangles (blue) denote points {Eini{N = 
114), M), (^inf(iV = 114), Me), and (^inf(iV = 114),^), re- 
spectively. Solid (red), dotted (green), and dashed (blue) 
lines are fitting curves for them, respectively. The triplet 
(Xi5 X2, Xs) denotes the tensor size. 



trust the M - E curve, and if we take an estimate of 
the energy E G (0.54, 0.55) upon which various previous 
works agree, we can conclude that M G (0.275,0.327). 
However, our result is clearly larger than previous es- 
timates: M = 0.205(15) from DMRG calculations^, 
M = 0.205(10) by GFQMC with SR calculations^^, and 
M — 0.19(2) by series expansioiP^. As shown in Fig. [oj 
the reason for this discrepancy may be the spatial inho- 
mogeneity caused by disentangling with small tensors. In 
DMRG calculations, to suppress the effect of the bound- 
ary condition as pinning the field on boundary sites, only 
on-site magnetization at the center of the system was 
usecP. To suppress the effect of incomplete disentangling, 
we also use on-site magnetization only at centers of ER, 
Mc. Figure [lO]plots the average on-site magnetizations at 
six centers of ER. If we assume the same condition for ex- 
trapolating M, we can conclude that G (0.232, 0.298). 
This value is significantly close to other estimations. It 
is probable that our present estimate of the magnetiza- 
tion is an overestimate owing to the intrinsic bias of the 
tensor network that generally favors states with less en- 
tanglement entropy. From this view point, estimating 
the magnetization at a position with larger ER may be 
more appropriate than simply taking the spatial average. 

However, the angle between magnetic moments of 
nearest neighbor sites converges, even when the tensor 



B. Spatially anisotropic triangular lattice of J2 < Ji 

We will report results of variational calculations for a 
spatially anisotropic triangular lattice. We only consider 
the case of J2 < Ji in this study. Our main interest lies 
in the robustness of the spiral magnetic ordered state. 



1. Tensor network 

As in the classical model, since the wave vector of a spi- 
ral magnetic ordered state may be incommensurate, we 
have to be careful about the periodicity in the variational 
wave function. Because of the finiteness of the unit cell 
in the tensor network, the wave vector of the magnetic 
ordered state is restricted. This restriction may cause a 
strong bias in variational calculations, in contrast to the 
case of an isotropic triangular lattice, where the unit cell 
of our ER tensor network is perfectly consistent with the 
three-sublattice ordered state. 

To weaken the finite-size effect of the unit cell, we in- 
creased the number of ER levels. In the MERA ten- 
sor network, the size of the unit cell increases exponen- 
tially by the number of ER levels. We have used a ten- 
sor network with two ER levels. The unit cell covers 
6 X 19 X 19 = 2166 sites. In addition, we make all tensors 
in the unit cell independent. Thus, the wave vector re- 
striction is relaxed compared to the tensor network with 
a single ER level. In detail, the reciprocal vector is writ- 
ten as 



/ 33 \ 

4332 V 63^3 J 



m 
4332 



41 
11^3 



(5) 



size is small. Figure 10 shows the average angle between 
magnetic moments of nearest neighbor sites, 0^ for vari- 
ous tensor sizes. The angle between magnetic moments 



where the number of independent sets, {l^m)^ is 2166, 
because of the skew periodic arrangement of the unit 
cell in Fig. [5] This number is 19 times that for the 
single ER case. Although the number of independent 
causal cones also becomes 19 times greater than before, 
we used parallel computing to calculate them. How- 
ever, because of the memory size limit in a computa- 
tional node, our calculation is limited to a tensor of size 
X = (Xi,X2,X3,X4,X5,X6) = (2,8,4,4,8,4). 
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2. Energy and quantum mutual information 

We have performed the variational calculations by 
PBC and infinite-size schemes with two ER levels {N = 
2166). 

As we see in the case of the single ER level, there is 
spatial inhomogeneity resulting from the structure bias 
of ER in both cases. By increasing the tensor size, we can 
systematically improve the quality of the tensor network 
states as before. In the single-ER calculation, as shown in 
Fig. [sj^d), a weak entanglement region between unit cells 
exists in the infinite-size scheme. However, it disappears 
in the infinite-size scheme with two ER levels. Because 
of the large unit cell obtained by the two ER levels, the 
finite-size effect of the unit cell for entanglement entropy 
is sufficiently removed. In the following, we will mainly 
report the results of an infinite lattice, which directly 
corresponds to the thermodynamic limit. 

Figure [TT] shows the energy per site for anisotropic 
cases from J2/J1 = 0.5 to 1.0. The tensor sizes in the 
tensor network with two ER levels are x = (2, 2, 2, 2, 2, 2), 
(2,4,4,4,4,4), and (2,8,4,4,8,4). The value of energy 
is improved by increasing the tensor size. In particu- 
lar, the results for the MERA tensor network are bet- 
ter than those of VMC calculations^^ in the region of 
J2/J1 > 0.75. Even at J2/J1 = 0.7, the result of 
MERA by using the infinite-size scheme with two ER lev- 
els (TV = 2166) is a little (0.5%) higher than that of VMC 
calculations^^. The difference may removed by the initial 
condition or by a tensor optimization process. However, 
the difference gets worse in the stronger anisotropic re- 
gion for J2/J1 < 0.7. Figure 12 shows the quantum mu- 



tual information (QMI) of nearest neighbor sites along 
the Ji and J2 axes. QMI represents the quantum corre- 
lation of two sites. If there is only a classical correlation 
between two sites, QMI is zero. The QMI of two sites i 
and j is defined as 



-0.44 
-0.46 
-0.48 
-0.50 
-0.52 
-0.54 
-0.56 



X 



VMC 

1 X = (2,2,2,2,2,2) 
% J =(2,4,4,4,4,4) 
# x = (2,8,4,4,8,4) 

4 X =(2,8,4,4,8,4) (PBC) 

m X =(2,8,8) 
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FIG. 11. (Color online) Energy per site by using tensor net- 
works with two ER levels (N — 2166). The results of the 
infinite-size scheme are mainly plotted. The results of the 
PBC scheme with two ER levels {N — 2166) and the infinite- 
size scheme with a single ER level [N — 114) are also plotted 
only for the largest tensor size. The VMC results are adapted 
from Tables I and III in Ref . ITOl 
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FIG. 12. (Color online) Quantum mutual information of near- 
est neighbor sites along Ji and J2 axes using the infinite-size 
scheme with two ER levels (N = 2166). The size of the ten- 
sors is X = (2,8,4,4,8,4). 



where Si and Sj are the entanglement entropy of site i 
and j, respectively, and Sij is the entanglement entropy 
of two sites i and j. As shown in Fig.[T2j when the spatial 
anisotropy increases (J2/J1 decreases), the QMI along 
the Ji and J2 axes increases and decreases, respectively. 
However, we assume an isotropic entanglement structure 
in our tensor network (see Fig. |3|. The mismatch may 
cause the poor performance in the stronger anisotropic 
region. Therefore, in the following, we will focus on the 
weak anisotropic region, J2/J1 > 0.7. 



3. Magnetization 

In the region 0.7 < J2/J1 < 1, the MERA tensor net- 
work states break SU(2) symmetry and they have finite 
on-site magnetizations for both finite and infinite-size lat- 
tices. Figure 13 shows magnetization and the average 



angle between magnetic moments on J2 links for three 
cases, J2/J1 = 0.7, 0.8, and 0.9, calculated using infinite- 
size scheme. As in the isotropic case, the dependence on 
tensor size remains. We cannot simply extrapolate them 
in the limit of the infinite dimension. However, even 
if the ground-state energy is about 0.1 lower than our 
best results, the extrapolated values from fitting curves 
are finite. Therefore, our results suggest that the ground 
states are magnetic. The wave function from the infinite- 
size scheme is a correct quantum state in the thermody- 
namic limit. Thus, at least, the magnetic state is a good 
candidate for the ground state in this model. In recent 
VMC calculations^^, the disappearance of magnetization 
for J2/J1 < 0.8 was reported. However, in our results, 
the magnetizations smoothly change even in the region 
of J2/ Ji < 0.8. 

The angle between magnetic moments of nearest neigh- 
bor sites weakly depends on the tensor size, as shown in 
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FIG. 13. (Color online) Average on-site magnetization and 
average angle between magnetic moments on J2 links. These 
results are obtained using the infinite-size scheme with two 
ER levels {N — 2166). The horizontal axis denotes en- 
ergy per site, Eini{N = 2166). The left and right verti- 
cal axes denote magnetization, M, and angle between mag- 
netic moments on J2 links, O2, respectively. Crosses, cir- 
cles, and stars denote (^inf(A^ = 2166), M) at J2/J1 = 0.9, 
0.8, and 0.7, respectively. Squares, triangle, and diamonds 
denote (^inf(A^ = 2166), 6*2) at J2/J1 = 0.9,0.8, and 0.7, 
respectively. The three curves are quadratic fittings for 
(Eini{N — 2166), M) points. The three lines are linear fittings 
for (Eini{N = 2166), ^2) points. For all cases of J2/J1, the 
tensor sizes of the left-most, middle, and right-most points are 
(Xi,X2,X3,X4,X5,X6) = (2,8,4,4,8,4), (2,4,4,4,4,4), and 
(2,2,2,2,2,2), respectively. 



Fig. 13 There are two groups of pairs of nearest neigh- 
bor sites: One is defined on the Ji links and the other 
is defined on the J2 links. We define 6i as the average 
angle between magnetic moments on Ji links. In detail, 
O2 may be split into two groups, 02a and ^26, which cor- 
respond to two directions along the J2 axis. First, in 
all results in the anisotropic region, the values of the 
sum Oi + 02a + 02b are 359.5(6). Thus, all magnetic mo- 
ments are coplanar, i.e., always lie on the same plane. 
Figure [M] shows the average angle O2 of MERA tensor 
network states. We plot results only for the largest ten- 
sor size. There is no discrepancy between 02a and 02b in 
the cases of two ER levels. As shown by the solid (red) 
circles and the solid (blue) triangles in Fig. [l4j the aver- 
age angle smoothly changes from 115.9(2) to 102(2) when 
J2/J1 decreases from 0.95 to 0.7. Thus the wave vectors 
of magnetic order are incommensurate. In contrast, as 
shown by the open (green) squares and the open (green) 
diamonds in Fig. [l4j the average angle suddenly changes 
around J2/J1 = 0.825 in the case of single ER. In detail, 
under J2/J1 < 0.8, O2 splits into 02a and 6>26. In other 
words, the reflection symmetry along the Ji axis breaks 
in J2/J1 < 0.8 by using the single ER level. Similar re- 
sults have been reported in a cylindrical lattice with a 
width as narrow as 6 in DMRG calculations^^. The rea- 
son for the behavior of angle in the case of single ER is 
the strong finite-size effect of the unit cell in tensor net- 
works. Thus, a large unit cell is necessary to capture the 



FIG. 14. (Color online) Average angles between magnetic 
moments on J2 links. The solid line denotes results of series 
expansion methods^^. The dotted line shows the value from 
the classical anisotropic triangular model. 



incommensurate state. The dashed line in Fig. [14] is the 
angle between magnetic moments on the J2 links for the 
classical spatially anisotropic antiferromagnetic Heisen- 
berg model. It is different from the results of MERA 
tensor networks. The solid line in Fig. [14] shows the re- 
sults of the series expansion method^^. Although the 
methods are very different, the series expansion results 
and those from the present work are in good agreement 
with each other. This fact gives a strong evidence for 
existence of a stable spiral phase. The change of wave 
number is larger than that of the classical one. Quantum 
fluctuations weaken the effective coupling between chains 
and enhance the incommensurability. 



IV. CONCLUSIONS 

Using ER tensor networks, we numerically studied the 
ground states of spin-| Heisenberg antiferromagnets on 
anisotropic triangular lattices. 

Since the area law of entanglement entropy holds in an 
ER tensor network, in principle, this new method may be 
effective for capturing high entanglement quantum states 
as expected in frustrated quantum magnets. Since we 
only assume an entanglement structure, our results will 
serve as a new piece of evidence independent of the pre- 
vious ones, with a totally new kind of bias. 

We reported numerical results using tensor networks 
with one and two ER levels, which correspond to = 114 
and N = 2166 unit cells, respectively. First, we con- 
firmed the 120° magnetic order ground state at the 
isotropic point Ji = J2 by using a tensor network with 
a single ER level {N = 114). The entanglement entropy 
was more sensitive to the direct product state on the top 
layer with a small unit cell size. Second, using the tensor 
network with two ER levels {N = 2166), we found a sta- 
ble spiral magnetic structure with incommensurate wave 
vectors at least in the anisotropic region 0.7 < J2/ Ji < ^^ 
In particular, the angle between magnetic moments on 
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nearest neighbor sites agrees very well with results of the 
series expansion method^^, which is a very different ap- 
proach. 

However, the spiral phase that we found overlaps the 
disordered phase reported in VMC calculations^^ Al- 
though we can roughly extrapolate magnetization in the 
limit of infinite dimension, we did not find the sharp de- 
crease in magnetization around J2/J1 = 0.85 reported in 
the VMC calculatioiffB. By increasing; the dimension of 
tensor indices, and by modifying the structure of the ER 
tensor network, we hope that we can obtain a complete 
answer in the near future. In particular, to overcome 
the computational cost, we may need to explore less de- 
manding methods in future studies such as combining the 
tensor network method with Monte Carlo sampling^^. 

In real materials such as /^-(BEDT-TTF)2 Cu2 (CN)3 
and EtMe3Sb[Pd(dmit)2]2, high-order interaction may 
play an import role. In particular, models with ring ex- 
change were discussed to explain the disordered behavior 



in real material^. The ER tensor network may be useful 
for studying ground states of such models. 
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